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EVALUATION 


1.  This  evaluation  covers  the  final  report  on  the  subject  contract 
by  Boston  College*  This  report  describes  the  theory  and  application 
of  several  methods  of  signal  analysis  and  data  processing  for  both 
ionospheric  HF  propagation  and  groundvave  LF  propagation. 

2.  The  objective  of  this  work  was  to  develop  advanced  signal  processing 
techniques  that  would  improve  the  target  detectability  of  an  AF  radar  in 
a  realistic  environment*  Several  methods  were  to  be  compared.  In  the 
low  frequency  part  of  the  spectrum  the  analysis  was  aimed  at  describing 
the  propagation  of  radio  waves  over  an  irregular,  finite  conducting 
ground. 

3.  The  results  of  the  radar  signal  processing  analysis  are  important 
although  they  are  relatively  negative  in  character.  That  is,  the 
advanced  methods  seem  to  produce  no  significant  improvement  over  the  more 
conventional  analysis.  This  conclusion  will  not  necessarily  be  true  of 

all  types  of  data.  The  effort  described  in  this  report  on  the  low  frequency 

2undwave  propagation  was  a  major  contribution  to  an  operational  tool  for 
an  type  systems. 
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1.  INTRODUCTION 


This  final  report  discusses  the  theory  and  application  of  various  types 
of  analysis  and  data  processing  for  both  ionospheric  HF  signal  propagation 
and  groundwave  LF  propagation. 

Much  of  the  analysis  pertaining  to  power  spectra  and  coherence  has  been 
dealt  with  in  previous  reports  [2].  The  results  of  applying  the  Maximum 
Likelihood  Method  to  aircraft  signal  data  has  been  detailed  in  [3] ,  [16] , [17] . 
The  main  effort  of  this  report  will  deal  with  the  Maximun  Entropy  Estimator 
Method,  its  derivation  and  application  to  data  from  the  Doppler  Arrival  Angle 
Spectral  Measurements  (DAASM)  program  and  to  the  sea  backscatter  measurements 
from  the  Ava- Dexter  HF  radar  facility  of  Rome  Air  Development  Center. 

Effort  was  also  directed  to  the  complimentary  task  of  decoding,  display, 
and  analysis  of  A.R. P.A.  Polar  Fox  II  ionospheric  backscatter  data,  represent¬ 
ing  Doppler  spectra  vs.  range. 

The  LORAN  grid  prediction  program  called  HUFIXJC  was  a  distinct  problem 

that  sought  to  combine  groimdwave  propagation  theory  with  a  suitable  model 

/ 

for  terrain  and  soil  electrical  characteristics  in  order  to  derive  an  adaptive 
estimate  of  time  of  arrival  differences  for  various  signals  used  in  the 
LORAN-D  navigation  system. 

An  extensive  effort  was  made  in  Section  5  to  extract  the  most  information 
regarding  radar  refraction  and  electron  density  profiles  from  a  survey  of 
vertical  incidence  ionograms,  using  existing  techniques  and  programs. 

2.  MAXIMUM  ENTROPY  ESTIMATOR 
2. 1  Analytical  Background 

This  study  is  intended  to  evaluate  the  performance  of  the  Fast  Fourier 
Transform  vs.  the  Maximum  Entropy  Method  estimate  in  the  processing  of  short 
data  sequences  in  an  OTH  radar  environment.'  Frequently,  a  radar  is  designed 
to  operate  using  spectra  estimated  from  samples  whose  duration  corresponds 
to  minimum  coherency  intervals.  Several  such  short  spectra  are  averaged 
before  decisions  are  made  concerning  the  presence  or  absence  of  target 
signals.  The  preference  for  the  frequency-domain  as  the  starting  point  of 
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analysis  is  dictated  by  the  presence  of  strong  groimd  clutter  returns  in  a 
frequency  band  around  DC.  The  target  signals  will  be  well  separated  from  this 
clutter  signal  because  of  the  Doppler  frequency  shift  due  to  the  notions  of  the 
target.  Since  the  question  of  determining  the  maximum  sample  size  for  evalua¬ 
tion  of  spectra  depends  crucially  on  the  coherence  times  and  scan-mode  of  the 
radar,  a  typical  specification  for  the  operational  mode  of  the  radar  and  the 
ionospheric  environment  which  effects  coherence  times  is  assumed. 

Once  the  time  samples  are  transformed  into  the  frequency  domain,  vectorial 
addition  of  corresponding  Doppler  Fourier  spectra  from  successive  blocks 
proves  futile  if  block  size  is  too  small.  Incoherent  averaging  or  scalar 
addition  of  corresponding  power  spectra  is  the  only  method  available  when  small 
block  sizes  (short  time  samples)  are  employed.  However,  if  higher  resolution 
spectra  are  calculated  from  the  same  temporal  sequences  appropriately  augmented 
with  zeroes,  the  spectra  can  then  be  incoherently  added  to  obtain  some  expected 
inprovement  in  signal  estimate  reliability  over  the  small  block  average  (a 
discrete  Fourier  Transform  (DFT)  that  computes  more  frequency  points  than  that 
of  the  time  series  is  equivalent  to  a  FFT  of  the  augmented  time  series) . 

Techniques  for  Signal  Processing  -  Two  techniques  for  signal  processing 
have  been  employed  on  data  using  actual  aircraft  signals.  They  are  respective¬ 
ly  incoherent  FFT  averaging  (Bartlett -Welch)  and  incoherent  maximum  entropy 
spectral  averaging.  Here  FFT  stands  for  Fast  Fourier  Transform.  The  different 
techniques  will  be  briefly  outlined  in  the  following  paragraphs. 

FFT  Incoherent  Averaging  -  In  the  incoherent  averaging,  a  data  sample 
of  say  128  complex  valued  points  is  broken  into  15  overlapping  segments  of 
16  points  each  and  a  FFT  is  performed  on  each  of  the  15  blocks  or  segments. 

An  estimate  of  the  power  spectra  is  made  for  each  of  the  blocks  and  an  average 
of  the  corresponding  Doppler  power  spectra  over  the  blocks  is  computed. 

Since  no  phase  information  is  utilized  in  computing  the  average  the  method 
is  called  incoherent  averaging. 

MEM  Incoherent  Averaging  -  The  second  technique  is  similar  to  the  first 
one  except  that  the  power  spectrum  for  each  block  is  calculated  by  using  the 
maximum  entropy  method  (M:M) .  The  MEM  is  well  known  and  has  been  described 


by  Burg  (4),  Lacoss  (1),  and  Smylie  et  al.  (5).  Varad  (2)  has  independently 
developed  and  applied  the  MEM  algorithm  to  evaluate  power,  coherence  and 
phase  spectra  where  the  input  samples  are  complex  valued  multi-channel  data 
points.  Superior  spectral  resolution  using  MEM  estimates,  particularly  for 
short  data  samples,  was  expected.  The  MEM  involves  recursive  estimation  of 
Anderson  coefficients  (5) .  The  MEM  spectra  are  averaged  in  the  same  manner  as 
the  incoherent  averaging  technique.  The  mathematical  details  for  these  tech¬ 
niques  are  presented  in  Appendix  A. 

2.2  Applications  to  DAASM  Data  [16] , [17] 

Aircraft  target  data  from  the  DAASM  program  presents  a  signal  displaced 
from  the  ground  clutter,  except  for  velocities  corresponding  to  the  Nyquist 
frequency,  ±4  Hz.  The  basic  Doppler  frequency  range  is  8  Hz.  A  data  set  is 
comprised  of  128  equally  spaced  complex  time  samples  from  each  antenna  channel 
and  each  range  gate.  No  interlacing  from  the  two  sets  of  equally  spaced  time 
samples  was  used  in  this  study,  so  that  Doppler  frequencies  greater  than  ±4  Hz 
are  aliased  into  the  basic  Doppler  frequency  range.  The  data  sample  is  divided 
into  15  blocks  of  16  points  each  with  a  50%  overlap  factor  to  reduce  variation 
from  block  to  block.  Figures  1  through  6  present  the  16  point  incoherent 
block  average  of  the  spectra  provided  by  the  MEM  with  16  and  256  frequency 
points  calculated,  compared  with  the  DFT  for  16  and  256  computed  frequency 
points.  The  MEM  was  calculated  for  2,  4,  and  7  filter  poles.  The  FFT  has  a 
3  point  Hanning  window  applied. 

2.3  Applications  to  Sea  Backscatter  Data 

Sea  scatter  data  from  RADC  Ava-Dexter  HF  radar  was  analyzed  with  both  the 
FFT  and  complex  MEM  techniques  using  a  variety  of  choices  for  data  length. 

When  reflected  from  a  land  mass,  the  spectra  typically  peak  at  zero 
frequency  (no  Doppler  shift),  indicating  no  local  movement  in  the  reflecting 
area.  However,  when  reflected  from  an  ocean  region,  the  returned  signals  yield 
two  distinct  Doppler  peaks  the  amplitude  of  which  is  proportional  to  the 
direction  of  local  ocean  wave  movement  in  the  reflecting  area.  Other  objects 
moving  through  the  reflecting  area,  such  as  ships  or  airplanes,  can  also  be 
detected  as  lower  level  peaks  within  or  near  the  same  Doppler  range  as  the  sea 
clutter. 
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Because  the  spectral  band  of  interest  extends  only  over  a  small  area  of 
the  lower  frequencies,  fine  resolution  is  required  in  the  frequency  domain 
to  readily  observe  the  two  Doppler  peaks  and  the  frequency  band  between  them. 
Typically,  very  long  time  samples  of  radar  backscatter  have  been  recorded  and 
processed  to  the  spectral  domain  by  conventional  Fourier  analysis.  Uhder 
these  conditions,  the  Doppler  peaks  are  resolved  directly  as  a  function  of  the 
length  of  the  time  samples  [15].  Theoretically,  the  same  spectral  data  resolu¬ 
tion  should  be  able  to  be  obtained  from  shorter  time  samples,  where  spectral 
analysis  would  be  performed  by  the  recently  developed  Complex  Maximum  Entropy 
Method. 

The  purpose  of  this  effort  then  is  to  verify  existence  of  certain  Doppler 
signals,  obtained  from  conventional  Fourier  analysis  of  long  time  data  samples 
read  from  magnetic  tape  recordings  of  ocean  radar  backscatter.  Having  shown 
the  resolving  power  of  this  method,  relative  to  data  length,  a  comparative 
spectral  resolution  analysis  by  Complex  MEM  should  theoretically  then  deter¬ 
mine  some  minimum  data  length  which  can  realistically  be  generally  applied 
to  signal  detection  from  radar  backscatter  information. 

Analytical/Programmable  Effort  -  Computer  programs  have  been  written  to 
read  given  magnetic  tape  recordings  of  time  domain  radar  backscatter  data. 

For  initial  observation  of  power  spectra  behavior,  computational  programs 
were  written,  and  different  data  sample  lengths  were  analyzed  by  use  of  Fast 
Fourier  Transforms  acting  on  the  raw  time  domain  data.  Initial  analysis  of 
the  power  spectra  showed  that  more  detailed  information  would  be  desirable. 

As  the  radar  information  extends  over  64  different  range  bins,  one  could 
choose  between  processing  by  the  so-called  Infinite  Fast  Fourier  Transform  the 
entire  64  range  samples  each  with  very  long  data  lengths,  or  stripping  informa¬ 
tion  by  range  and  then  processing  time  sequences  for  the  selected  ranges.  The 
former  required  large  computer  core  and  processing  time,  thereby  reducing 
computer  turnaround  time;  the  latter  required  writing  additional  data  handling 
programs  with  intermediate  disk  storage  of  information.  As  it  was  felt  that 
the  Complex  MEM  analysis  would  probably  not  be  necessary  over  all  64  ranges,  the 
second  option  was  undertaken. 
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A  program  was  written  from  which  disk  files  of  time  domain  information 
were  generated  to  allow  selection  of  specific  ranges  by  the  main  analytical 
routine.  For  increased  spectral  resolution,  an  FFT  was  applied  to  the  zero- 
augmented  time  series  with  a  cosine  squared  window  (Hanning  filter)  imposed 
only  on  the  basic  input  time  interval.  [Bart let t- Welch  block  averaging  was 

applied  to  sequences  of  computed  power  spectra  to  improve  overall  signal 
estimation  reliability.]  The  resulting  power  spectra  are  displayed  in  Figures 
7  and  8. 

2.4  Results  6  Conclusions 

The  new  MEM  routine  produced  reasonable  and  consistent  Bragg  line  spectra 
(Figs.  9-12).  Power  spectra  from  the  FFT  were  almost  identical,  indicating 
that  the  natural  line  width  was  being  displayed  at  least  for  a  time  series  of 
128  points  or  more. 

The  anticipated  order  of  magnitude  of  Signal -to-Noise  improvement  was 
not  observed  when  comparing  properly  windowed  Fourier  spectra  to  optimal  MEM 
spectra.  Here,  optimal  is  used  in  the  sense  that  some  optimum  number  of 
filter  poles  was  utilized  to  produce  a  spectrum  of  minimum  noise  and  minimum 
smoothing.  While  some  MEM  spectra,  produced  with  larger  numbers  of  filter 
poles,  appear  to  have  more  sharpness  in  spectral  peaks,  the  same  sharpness 
is  apparent  in  the  noise  level.  Overall  frequency  resolution  was  comparable 
when  comparing  spectra  generated  from  input  of  the  same  data  lengths.  MEM 
analysis  of  shorter  data  lengths  did  not  produce  significantly  observable 
spectral  peaks  in  the  target  areas.  Overall,  target  detectability,  small 
signal  detection  in  the  higher  Doppler  frequency  ranges,  was  not  very  much 
different  between  the  Fourier  and  MEM  spectra. 

One  conclusion  to  be  drawn  from  this  investigation  is  that  it  seems 
certain  classes  of  data  are  not  of  a  nature  as  to  allow  spectral  analysis 
methods  based  upon  prediction-error  filters  to  fully  optimize  data-adaptive 
features.  The  nature  of  such  information  is  not  defined  within  the  context  of 
this  study.  A  possible  area  of  investigation  might  be  system  bandwidth 
characteristics. 

Similarly  for  the  DAASM  aircraft  data,  the  MEM  with  a  256  point  DFT 
offers  only  slight  improvement  in  resolution  over  the  equivalent  256  point  FFT 
of  the  zero- augmented  16  point  time  series;  compare  Figures  4  and  6  or  Figures 


2  and  3.  In  general ,  the  number  of  spectral  peaks  produced  equals  the  number 
of  MEM  filter  poles. 


3.  HF  IONOSPHERIC  PROPAGATION 
3. 1  Data  Processing  S  Mode  Identification 

The  Polar  Fox  II  program  of  the  Defense  Advanced  Research  Projects  Agency 
(ARPA)  studied  ionospheric  backscatter  signals  over  the  6  to  26  Mlz  frequency 
range  in  the  northern  latitudes,  with  the  receiver  and  transmitter  located  at 
Caribou,  Maine  and  scanning  azimuths  from  -30°T  to  +60°T  during  1972.  The 
magnetic  tapes  generated  contain  two  types  of  data;  a  fixed  frequency,  azimuth 
scanning  mode  and  a  fixed  azimuth,  frequency  scanning  mode.  The  frequency 
scanning  mode  produced  ionograms  (frequency  vs.  range)  at  3  fixed  azimuths, 
while  the  azimuth  scanning  mode  gave  a  32  point  Doppler  spectra  at  230  range 
bins  (from  624  km  to  4074  km)  at  3  discrete  frequencies.  Figure  13  shows  the 
range  coverage  and  geography  of  the  experiment  and  Figure  14a  presents  the 
digicoder  display  generated  from  unpacking  and  decoding  the  frequency  scanned 
ionogram  data.  Reports  from  the  Raytheon  Co.  and  Lincoln  Laboratories  were 
the  sources  for  the  program  developed. 

The  fixed  frequency  backscatter  (B/S)  data  was  decoded,  unpacked  and 
calibrated  from  the  appropriate  hourly  calibration  signal;  a  noise  level 
determination  is  available  for  every  minute  (i.e.,  every  azimuth)  of  the  34 
minute  data  interval  and  this  is  used  as  a  threshold  for  displaying  the 
return  Doppler  signal.  Four  types  of  displays  were  developed: 

1.  Received  power  calibrated  in  dBw  from  the  reference  level  and 
thresholded  at  6  dB  above  the  4th  decile  of  noise,  32  Doppler 
bins  per  every  3rd  range  bin,  one  display  for  each  azimuth  beam 
position.  (Figure  14b) 

2.  A  "contour”  plot  of  this  display  with  a  different  letter  assigned 
to  every  10  dB  step. 

3.  A  listing  line  plot  of  amplitude  vs.  range  bin  for  the  0,  ±8 
and  16  Hz  Doppler  lines. 

4.  A  three  dimensional  plot  of  Doppler  vs.  range  vs.  power  amplitude 
for  any  selected  azimuth  or  frequency.  (Figure  15). 
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Together  these  plots  facilitate  the  propagation  mode  identification  with 
the  presence  of  ground  clutter,  interference  lines,  etc.  The  four  transmitter 
powers  recorded  indicate  both  total  power  and  antenna  operating  mode;  the 
azimuth  beam  pattern  gain  and  beamwidth  are  derived  from  the  published  system 
description  [19]. 

3.2  Backscatter  Reflectivity  Model 

The  next  step  is  to  apply  the  two  basic  calculations  of  volumetric 
backscatter  (VBR)  and  surface  backscatter  reflectivity  (SBR)  to  this  total 
received  power.  The  takeoff  elevation  angle  is  determined  as  a  geometric  func¬ 
tion  of  range;  ionospheric  absorption  as  a  function  of  geomagnetic  latitude 
[7],  and  the  assumed  elevation  orthogonality  angle  are  also  needed  (Appendix 
B) .  The  VBR  and  SBR  equations  were  developed  by  G.  S.  Sales  [6].  The 
resulting  calculations  can  be  made  for  ionospheric  E  and  F  layer  heights  and 
different  propagation  mode  options.  (Figures  16  and  17,) 

The  synoptic  data  tapes  were  revived  and  used  to  display  larger  periods 
of  mode  identifications;  a  short  program  was  coded  for  this  purpose.  Previous 
attempts  to  analyze  this  data  at  other  facilities  has  produced  results 
questionable  from  the  point  of  view  of  mode  identification  and  which,  moreover 
took  no  account  of  antenna  sidelobes  in  the  azimuth  identification.  The 
continuing  analysis  effort  will  attempt  to  incorporate  antenna  radiation 
patterns,  at  least  to  the  extent  of  avoiding  gross  mislocation  of  the  clutter 
scattering  area.  This  is  a  definite  possibility  in  view  of  the  poor  sidelobe 
suppression  as  the  antenna  is  steered  to  large  angles  off  boresight.  An 
effort  must  be  made  to  check  the  internal  consistency  of  the  data  for  clues 
to  sidelobe  or  backlobe  reception  as  well  as  reasonable  agreement  with  nominal 
reflectivity  values  for  standard  types  of  clutter.  Also,  samples  of  data 
recorded  at  times  and  frequencies  differing  by  a  small  amount,  showing  conflict¬ 
ing  mode  identification  for  the  same  geographical  region,  will  be  studied  more 
closely.  Clearly,  error  can  result  if  the  received  backscatter  came  from  a 
different  direction  than  that  immediately  indicated  by  the  antenna  bearing, 
since  certain  assumptions  must  be  made  about  the  ionosphere  in  the  region  from 
which  the  backscatter  is  assumed  to  originate. 
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4.  LF  GROIMDWAVE  PROPAGATION 


4.1  System  Description  -  LORAN-D 

LORAN-D  is  a  hyperbolic  system  of  navigation  wherein  a  radio  navigation 
fix  is  obtained  as  the  intersection  of  two  hyperbolic  lines  of  position 
(Fig.  18).  These  hyperbolic  lines  are  synthesized  from  a  set  of  three  pulse 
transmitters  comprising  a  master  M  and  two  slaves,  S ^  and  S2  separated  by 
a  convenient  baseline  length  d^  and  d^-  The  navigator's  line  of  position  is 
expressed  as  a  constant  time  difference  TD  where 

nl 

TD  *  —  (d.  +d  -d  )  +  t  d.  ♦  t  d  -  td  ♦  e„ 
c  o  s  nr  cT>  c  s  cm  s 

where  n^  =  index  of  refraction  of  air 

c  =  velocity  of  light 

d^  =  length  of  baseline,  geodetic  line  from  master  to  slave 

dm  »  geodetic  line  from  observation  point  to  master  transmitter 

dg  =  geodetic  line  from  an  observation  point  to  a  slave  transmitter 

es  *  coding  or  emission  delay  (retransmission  time  of  the  slave  signal) 

t  »  phase  correction  for  propagation  over  a  path  of  length  d^dm. 
c  dg,  secs. 


The  distances  dg,  dm  and  d^  are  also  geodesics  to  a  particular  point  0 
along  AB  or  between  master  M  and  slaves  s^  and  s^.  Curved  lines  have  been 
used  to  note  that  geodesics  are  characterized  by  a  continuous  azimuth  or 
direction  change  along  the  line.  Two  hyperbolic  lines  of  position  TD1  and  TD2 
determine  the  point  0  at  their  intersections.  Given  the  geographic  coordinates 
of  M,  Sj  and  S the  coordinates  of  0  can  be  calculated.  However,  the  phase 
correction  expressed  in  seconds,  t  ,  must  be  known  for  precision  determination. 
The  phase  correction  becomes 

WtC 

<b  =  — r  radians 


f  =  w/2tt  Hertz 


This  phase  is  in  turn  composed  of  a  primary  delay  due  to  a  homogenous 
terrain  and  a  secondary  delay  due  to  irregularities  and  inhomo geni ties  in  the 
terrain.  The  primary  phase  factor  along  the  geodetic  distance  is  calculated 
from  the  Sodano  method  [11].  The  secondary  phase  factor,  in  the  method  of 
Johler  and  Hyovalti  is  calculated  from  assumed  electrical  characteristics  of 
the  terrain.  Sodano' s  method  has  been  adapted  to  incremental  geodesic  steps 
by  Huf  ford  [8] ,  [9] ,  [10] . 

In  the  final  software  package,  the  overall  program  called  HUFLOC  accepts 
topographical  and  geological  information  from  the  random  access  data  base; 
applies  impedances  calculated  from  this  data  base  to  the  geodesic  signal  paths 
determined  from  transmitter  and  receiver  coordinates,  calculates  the  time  of 
arrival  of  each  groundwave  (combining  primary  and  secondary  phase  factor),  and 
finally  calculates  the  two  required  time  differences  TD1  and  TD2.  The  two 
basic  programs  include  a  means  of  calculating  geodesic  distances  between  trans¬ 
mitter  and  receiver  by  Huf  ford's  method,  and  then  proceed  to  perform  the 
integration  of  the  signal  propagation  over  the  geodesic  distances  to  yield  the 
components  (amplitude  and  phase)  of  the  waveform.  A  combined  integral  equation 
and  classical  electromagnetic  technique  is  used  [8],  [9],  coded  in  subroutine 
INEQ2E. 

4.2  Terrain  Electrical  Characteristics 

This  effort  involved  the  establishment  of  the  random  access  data  base 
working  from  a  set  of  terrain  elevation,  geologic  age  and  soil  classification 
tapes  supplied  through  the  Initiator,  the  LORAN  S.P.O.  office.  A  total  input 
of  37  tapes  is  used  to  generate  a  data  base  grid  of  691,200  points  containing 
packed  values  of  elevation  and  either  complex  impedance  or  geology  and  soil 
code,  located  every  30  arc-seconds  over  a  6°  by  8°  geographic  area. 

Each  input  geology  code  tape  was  itself  generated  in  three  steps: 

1.  Conversion  from  digitized  table  coordinates  to  geographic  coordinates, 
corrected  for  longitude  line  conversion 

2.  A  SORT- MERGE  routine  to  select  and  sequence  the  resulting  data 

3.  A  final  conversion  to  the  interpolated  30  arc-second  grid. 

This  conversion  sequence  was  developed  and  coded  for  use  here  on  the  CDC  6600. 
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The  actual  coaplex  i ape dance  is  derived  froa  the  correspondence  table 
of  actual  conductivities  and  dielectric  constants  along  with  an  assuaed  aodel 
of  the  soil  and  rock  strata;  the  resulting  inpedance  is  overwritten  on  the 
fundamental  data  base  as  altitude  and  phase.  Ten  bits  are  allocated  to 
elevation,  9  bits  to  aaplitude  and  11  bits  to  phase  produced  an  aaplitude 
accuracy  of  ±2  ohas,  phase  accuracy  of  ±.001  radians  and  elevation  to  ±7 
aeters.  These  allocations  were  chosen  on  the  basis  of  the  sensitivity  of  the 
coaputed  tiae  of  arrival  of  the  groundwave  and  the  aeasureaent  accuracy  of 
the  conductivity  aaps.  The  soil  strata  aodel  chosen  can  be  changed  either 
for  geophysical  reasons  or  as  a  result  of  the  sensitivity  and  estiaation 
procedures  outlined  below. 

4.3  Time  Difference  Grid  Prediction 


The  combined  HUFLOC- INEQ2E  program  accesses  the  ran don  access  data  base 
by  means  of  an  index  generated  froa  input  path  vector  increaents  along  the 
geodesic  from  Master-,  Slave  1-  and  Slave  2-to-receiver  and  proceeds  to  solve 
the  integral  equation  for  field  intensity  and  phase  delay,  primary  and 
secondary.  A  complete  description  of  the  prograa  is  given  in  [12]. 

Next,  a  least-squares  estimate  with  a  Kalman  filter  was  applied  for  the 
empirical  correction  of  the  soil  conductivity  aodel.  The  specific  approach 
developed  was: 

a.  The  two  values  for  topsoil  conductivity  that  occurred  most 
frequently  in  the  geographical  area  of  interest  would  be  identified. 

b.  The  random  access  data  base  would  be  created,  and  the  program  which 
computed  the  transmission  time  delay  would  be  run  to  produce  cal¬ 
culated  values  for  approximately  SO  locations  for  which  measured, 
airborne  data  were  available. 

c.  The  variation  of  transmission  time  delay  due  to  variation  in  topsoil 
conductivity  would  be  computed  by  changing,  in  turn,  each  of  the 
two  conductivity  values  where  they  occurred  in  the  geographical 
area  of  interest  and  by  proceeding  with  the  steps  outlined  in  para¬ 
graph  b .  above . 

d.  The  differences  between  the  observed  and  calculated  data,  the  rates 
of  change  in  transmission  time  delay  due  to  changes  in  each  con¬ 
ductivity  value,  estimates  of  the  error  in  the  initially  chosen 
values  of  conductivity,  and  estimates  of  the  error  in  the  observa¬ 
tions  would  be  input  into  the  Kalman  filter  program  to  produce  an 
improved  estimate  of  the  values  of  conductivity. 
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e.  Using  the  improved  values  of  conductivity,  the  steps  outlined  in 
paragraph  b.  above  would  be  repeated,  and  the  new  calculated  values 
could  be  contrasted  with  the  earlier  transmission  time  delay  values. 

In  order  to  select  the  two  most  frequently  occurring  values  of  topsoil 
conductivity,  a  modified  version  of  the  data  base  creation  program  was  written. 
This  program  displayed  all  the  values  of  topsoil  conductivity  in  a  map-like 
fashion.  It  was  immediately  obvious  that  there  were  extraneous  values  of 
conductivity  along  the  same  line  of  longitude  in  numerous  places  in  the 
geographical  area.  Many  of  these  values  were  corrected  by  scanning  along  the 
same  latitude  line  and  changing  those  values  which  differed  from  that  value 
wtich  occurred  immediately  before  and  after.  A  similar  check  was  perfoxmed 
for  the  values  of  bedrock  conductivity.  After  this  preliminary  work,  the  two 
most  frequently  occurring  values  of  topsoil  conductivity  were  identified. 

The  locations  in  the  random  access  data  base,  where  these  conductivity  values 
occurred,  were  identified  in  order  to  reduce  by  70%  to  90%  the  time  required 
to  cnange  the  data  base. 


Calculated  values  of  transmission  time  delay  were  made  for  43  locations 
for  which  observed  data  was  also  available.  Numerical  approximations  to  the 
partial  derivatives  of  transmission  time  delay  with  respect  to  each  value  of 
conductivity  were  made  by  varying  each  value  of  conductivity  ♦10%  and  -10%. 

New  data  bases  were  produced,  and  calculated  values  of  transmission  time  delay 
were  made.  The  partial  derivatives  were  approximated  by  the  quotient  of  the 
difference  in  transmission  time  delays  divided  by  the  difference  in  conduc¬ 
tivity  values.  These  calculations  involved  160  lengthy  computer  runs. 


The  variance  associated  with  each  data  measurement  was 

2 


o  2 
o  sinor 


2  . 


where  is  an  input  quantity  and  a  is  an  angle  formed  by  the  three  points., 
master  transmitter,  receiver,  and  slave  transmitter.  The  standard  deviation 
oq  of  the  measurements  is  assumed  to  be  approximately  100  or  200  nanoseconds. 
The  filter  program  was  run  several  times  in  order  to  determine  what  changes 
resulted  in  the  estimate  of  the  values  of  conductivity.  The  initial  values 
foT  the  state  vector  of  conductivity  values  and  the  state  covariance  matrix 
were  respectively: 


(02061 A  / 1.  x  10-6  0  \ 

008439 )  \  0  1.x  10 -6I  . 

Using  values  of  oq  equal  to  200  nanoseconds  and  100  nanoseconds,  respectively, 
the  following  estimates  of  the  state  vector  resulted: 


New  data  bases  were  produced  for  each  of  the  sets  of  values,  and  trans¬ 
mission  time  delays  were  calculated  for  six  of  the  40  locations  available.  In 
general,  the  differences  between  observed  and  calculated  time  delays  decreased 
as  compared  to  the  differences  obtained  using  the  original  data  base.  The 
second  estimate,  as  anticipated,  resulted  in  the  smallest  differences  of  the 
original,  first  changed,  and  second  changed  data  bases. 

4.4  Results 

The  predicted  time  difference  grid  in  general  was  about  200  nanoseconds, 
(equivalent  to  a  C.E.P.  of  about  60  meters)  different  from  measured  recon¬ 
naissance  values  from  aircraft  reception  at  about  1000  ft.  altitude  and 
somewhat  better  on  ground  reception:  the  core  memory  required  to  run  the 
HUFLOC  program  is  about  137K  and  the  execution  time  for  a  typical  target  is 
130  secs. 

5.  RADAR  REFRACTION  USING  VERTICAL  INCIDENCE  I0N0GRAMS 

The  data  from  the  Digisonde-128  were  stored  on  1200  reels  of  tape. 
Originally,  one  day's  worth  of  data  was  put  on  one  tape;  later  on,  two  days 
worth.  The  data  themselves  represent  a  digital  ionogram  which  covers  a 
frequency  range  from  0  Mlz  to  about  12  Wz,  although  this  upper  limit  can  be 
increased  to  16  KHz  or  even  higher.  On  the  vertical  scale  the  ionogram  has 
128  range  bins,  the  first  bin  corresponding  to  60  km.;  each  successive  bin 
represents  1.5  km.  On  the  horizontal  scale  the  step  size  is  generally  ,05 
Mlz  for  data  gathered  at  the  beginning  of  an  hour;  and  either  .10  or  .025  Hiz 
otherwise.  If  we  take  the  128  range  bins  and  multiply  times  the  12/. 05 
(*240)  frequency  positions  we  get  about  30,000  points  (or  possible  echoes) 
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for  each  ionogram.  Since  we  have  about  9  ionograms/hour,  and  the  Digisonde 
works  24  hour/ day,  it  is  apparent  that  the  amount  of  data  is  considerable. 
However,  it  turns  out  that  only  about  200  of  these  points  actually  represent 
the  traces  of  the  ionogram;  everything  else  is  "noise",  or  unwanted  signals. 
Therefore,  by  filtering  we  were  able  to  compress  the  data  from  1200  tapes 
down  to  about  10.  The  actual  method  of  filtering  is  described  in  [13]. 

The  compression  operation  also  involves  packing  the  data  for  each  echo 
into  half  of  one  CDC  word.  The  data  stored  include  the  number  of  the  fre¬ 
quency,  the  number  of  the  range  bin,  the  amplitude  of  the  echo,  and  the 
total  number  of  echoes  at  the  given  frequency.  This  information  requires 
28  bits  altogether.  Thus ,  the  data  for  two  echoes  can  be  stored  in  me  CDC 
word. 


Trace  identification  involves  the  determination  of  certain  critical 
frequencies  and  heights  from  the  ionograms.  The  actual  description  of  the 
techniques  involved  can  be  found  in  [13]. 

An  important  application  of  filtered  ionogram  data  is  that  of  construe  - 
ing  an  electron  density  profile.  The  programs  already  exist  for  this;  all 
it  is  necessary  to  do  is  to  select  a  group  of  points  from  the  ionogram  and 
feed  them  into  a  profile  program  and  get  a  plot  of  height  versus  electron 
density.  The  only  problem  is  selecting  the  points  properly.  This  is 
discussed  in  [14]  by  Cormier  and  Dieter.  Once  an  election  density  profile 
for  an  ionogram  is  calculated,  it  is  of  some  interest  to  apply  the  curve  in 
computing  refraction  corrections  for  a  radar  system.  The  problem  centered 
around  using  different  curve- fitting  routines  for  different  parts  of  the 
profile.  We  employed  a  spline  curve  over  the  lower  part  of  the  profile  and 
an  exponential  curve  over  the  upper  part.  This  latter  curve-fit  was 
described  in  a  previous  Boston  College  report  by  Power  and  Mart in e  [18]. 
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APPENDIX  A 

Maxima  Entropy  Method 


In  this  section,  the  Mathematical  derivations  will  be  presented  for  the 
different  techniques. 


1)  Incoherent  FFT 

Let  Xj(t),  n^Ct) . .  .x^t)  be  the  N  complex  time  samples,  e.g.,  N=128. 
The  N  sample  points  are  divided  into  L  blocks  of  M  sample  points  and 
a  FFT  of  length  M  is  performed  for  each  of  the  L  blocks.  Thus,  the 
FFT  is  given  by 

XK(w)  *  V  xn(t)  e_jwkn  , 

n*o 

where  K  is  the  discrete  frequency.  The  average  power  in  the  ktk 
frequency  is  calculated  from 

l  L  2 

sk(w)  *  r  E  • 


A  Hanning  filter  could  be  implemented  by  replacing  Xg(w)  with 
y  f v('*)  ,  i  y(^)  ,  i  v(^0 

V  J  '  4  Vl  2  *K  4  K+l 

The  end  points  X^(w)  and  X£j_j  are  wrapped  around  for  this  estimate. 


2)  Incoherent  MEM 

The  maximum  entropy  method  of  spectral  evaluation  seeks  to  determine 
coefficients  of  a  filter  whose  output  would  be  ideally  white  noise 
when  the  input  is  the  process  under  consideration.  Thus,  if  S(w)  is 
the  power  spectrum  of  the  input  data  and  H(w)  is  the  transfer  function 
of  the  filter,  then 

S(w)|H(w)|2  -  o2 

2 

where  a  is  the  variance  of  the  output.  The  coefficients  comprising 
H(w)  for  any  given  order  of  the  filter  are  chosen  so  as  to  minimize 
o2.  For  a  (K*l)tk  order  filter 


H(w)  ■  1  ♦  j  e"^w  *  ♦  a^  j  e'^w*^  ♦  .  ..a^  K  e"^w^ 


1  ♦  2  a„  .  © 
k-1  *»* 


so  that 


|1  ♦  2  aKk**H 
k=  1  K,K 

The  coefficients  a„  ,  are  calculated  in  a  recursive  manner  beginning 

*»K  2  2  2  2 
with  a11,a22...aK  j  K  j  and  aR  likewise  from  oQ  ,ak  . . .oK  Initially 


is  calculated  from 


V  =  P1 


N-l 

P1  "  1  Xi*i 

i»0 


th 


where  xj  represents  the  complex  conjugate  of  xk»  the  i  time  sample 


Next,  the  autocorrelation  for  zero  lag  rQ  is  set  equal  to 

ro  =  P1 


Now  r^,  the  autocorrelation  for  unit  lag  is  calculated  from  the  matrix 
equation 


where 


ro  *1  1 


rl  ro  all 


P2  "  °1  "  .M>*i  *  allXi+l^  *  I* 

1*0  l 


i  *  *ii  j 


and  is  obtained  by  minimizing  P2 


=  0 

3all 

Here  P2  is  being  evaluated  by  forward  and  backward  averaging  and  the  corres¬ 
ponding  ajj  is  given  by 

an  ■  -2  '  E(xixi  *  xi*i  xi‘i>) 

and  P2  =  (1  -  a^Jj) 

The  next  recursion  is  written  as 


xo x2  rn 

n  ro  rl  *  a22<l*h 

x2  t;  ro  (?)  I1 


0  ‘»21<o 


P3  is  found  in  a  similar  manner  as  P2  by  minimizing  with  respect  to  a22 
and  given  by 

P3  (1  a22  a22)  P2 


*21  *  all  +  a22  all 


N-3  ,  v 

-2  .fo  |Xi+1  *  auxi+2  Xi+2  -  au  Xi+1J 

{K*‘  *  *UXi*^  *  ^Xi*2  *  *nXi.J  } 


Further  recursions  can  be  similarly  performed  and  the  coefficients  at 


the  step. 


aKl,aK21*‘ 


aK,K-l,aK,K 

can  be  calculated  and  inserted  in  the  equation  for  the  power  spectra,  and 
of  course 


a 


2 

K 


P 


K*1 


is  also  available  recursively. 
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APPENDIX  B 


Three  basic  propagation  modes  are  considered,  two  associated  with 
ionospheric  backscatter  and  the  third  with  ground  backscatter  (Fig.  16) . 
Permutations  of  two  basic  heights,  100  km  and  300  km,  are  used  to  establish 
the  possible  propagation  modes  and  scattering  heights.  The  atmospheric 
absorption  loss  term  L  was  calculated  from  the  method  developed  in  [7]. 

The  ground  surface  reflectivity  pQ  is  defined  in  terms  of  measured 
quantities.  For  volume  scattering,  the  concept  of  orthogonal  scattering  with 
cylindrically  shaped  irregularities  has  been  introduced.  As  in  the  case  for 
a  surface  reflector,  the  azimuthal  beam  width  and  pulse  width  define  two  of 
the  three  required  dimensions  for  scattering  volume  in  the  ionosphere  which 
contributes  to  the  received  power  at  a  specific  time  dealy  or  range,  R.  The 
third  dimension,  R(li),  has  been  determined  using  ray-trace  techniques 
in  F-region  models  of  the  ionosphere  and  the  required  orthogonality  with  the 
earth* s  magnetic  field  for  backscatter. 

The  effective  vertical  angle,  ^H),  is  shown  in  Figure  17  and  is  deter¬ 
mined  by  ray  tracing  as  having  a  nominal  value  of  0.5  degrees. 


The  two  intrinsic  reflectivities  derived  by  Sales  [6]  to  characterize 
the  scattering  medium  are  the  surface  reflectivity 


3  4  2  ^ 

.  (4 *r  k  (l  ;  p 

ptgtgra2(area)  r 


units  o 


f  M2 /M2 


and  volume  reflectivity 


3  4  2  2N+1 

(**)  k  (L  j  p 

PtGtGrA2 (VOLUME)  R 


units  o 


f  M2/^ 


where  the  illuminated  area  =  t(R$)  in  M2  and  the  illuminated  volume  =  x(R4>) 
IR@]  in  and 
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T  =  pulse  width 

$  ■  azimuthal  beamwidth 

=  effective  vertical  angle 

N  =  0,  one  leg  mode 
■  1,  three  leg  mode 
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Figure  1.  DAASM  142-2,  2  Pole  MsM  vs.  FFT 
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Figure  5.  DAASM  142-4,  7  Pole  MEM  vs.  FFT 
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Figure  7.  Sea  Scatter,  128  Pt.  Time  Series  (Hanning 
Window)  512  Frequency,  FFT 
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Figure  8.  Sea  Scatter,  128  Pt.  Time  Series  (No 
Hanning)  S12  Frequency  FFT 
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AZIMUTH  =  40°  T 


Figure  14a.  Signal  Frequency  vs.  Range  Ionogram 
b.  Doppler  Frequency  vs.  Range  B/S 


72/0325/18  UT 
16.1  MHZ 
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GROUND  REFLECTION  AND  IONOSPHERE  1C  BACKSCATTER 

(THREE  LEG) 


Figure  17.  Effective  Vertical  Angle 
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